
quietly {
global pops=1
global impacts1=1
global impacts2=1
global impacts3=1

clear
set obs 1
capture drop wnb*
	forvalues y=1/$pops {
	forvalues z=1/${impacts`y'} {
	gen wnb`y'_`z'=.
	format wnb`y'_`z' %15.0f
	} 
	gen wnb`y'=.
	format wnb`y' %15.0f
	}
	gen wnb=.
	format wnb %15.0f

forvalues i=1/500 { //number of MC simulations
	if mod(`i',5)==0 noisily	di "`i'"
	replace wnb=0 in `i'

	*General weighting parameters
	rtriangle, min(54998) mode(61109) max(67219)
	global med = $x
	rtriangle, min(1.2) mode(1.6) max(2)
	global epsinc = $x
	global bnmin = 10
	global bnmax = 10
	global thresh_lo=.5
	global thresh_hi=5
	rtriangle, min(1.555) mode(4.212) max(6.873)
	global alpha=$x
	rtriangle, min(.067) mode(.131) max(.195)
	global beta=$x

	*Population 1 parameters
	global hhdist="& shp_hh==1"
	global unif=0 
	rtriangle, min(9000) mode(10000) max(11000)
	global min1 = $x
	rtriangle, min(180000) mode(200000) max(220000)
	global max1 = $x
	rtriangle, min(45000) mode(50000) max(55000)
	global modepop1 = $x
	
	*Impact 1_1 parameters VSL
	global pp_units1_1=1	
	global cost1_1=0
	*rtriangle, min(22500) mode(30000) max(37500)
	global modepp1_1 = $min1
	*rtriangle, min(7500) mode(10000) max(12500)
	global minpppct1_1 = 1 
	rtriangle, min(.25) mode(.5) max(.75)
	global maxpppct1_1 = $x
	rtriangle, min(3589000) mode(10875000) max(18169000)
	global avgwtp1_1 = $x
	rtriangle, min(.287) mode(.519) max(.75)
	global epswtp1_1 = $x
	rtriangle, min(63432) mode(70480) max(77528)
	global avgwtpinc1_1 = $x
	global avgwtp_simpop1_1 = 10875000 //does not vary, used for simpop
	rtriangle, min(9787500) mode(13050000) max(16312500)
	global nb1_1 = 13050000

	*compute wnbi for all impacts for all populations
	global temp=0
	forvalues y=1/$pops {
		global wnbtot=0
		global temp=${modepop`y'}
		global modepop=${modepop`y'}
		global min=${min`y'}
		global max=${max`y'}
		*global unif=${unif`y'}
		global hhdist="${hhdist`y'}"
		forvalues w=1/${impacts`y'} {
		global cost=${cost`y'_`w'}
		global modepp = ${modepp`y'_`w'}
		global minpppct = ${minpppct`y'_`w'}
		global maxpppct = ${maxpppct`y'_`w'}
		global avgwtp = ${avgwtp`y'_`w'}
		global avgwtp_simpop = ${avgwtp_simpop`y'_`w'}
		global epswtp = ${epswtp`y'_`w'}
		global avgwtpinc = ${avgwtpinc`y'_`w'}
		global nb=${nb`y'_`w'}
		wnb ${pp_units`y'_`w'}
*	} // temporary, for doing non-MC sensitivity analysis
*	} // temporary, for doing non-MC sensitivity analysis
		if $cost==0 {
			quietly replace wnb`y'_`w'=$wnb in `i'
			global wnbtot=$wnbtot+$wnb
			}
		if $cost==1 {
			quietly replace wnb`y'_`w'=$wnb*-1 in `i'
			global wnbtot=$wnbtot-$wnb
			}
		quietly replace wnb`y'=$wnbtot in `i'
	} // close w
	quietly replace wnb=wnb+wnb`y' in `i'

	} // close y
local temp=`i'+1
set obs `temp'
} // close i
*drop if missing(wnb1_1-wnb)==1
format wnb %12.0f
quietly sum wnb
noisily di "wnb=`r(mean)'"
drop if wnb==.
} // close quietly
di "Varying weighting parameters."